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Abstract. The TASEP is a paradigmatic model from non-equilibrium statistical 
physics, which describes particles hopping along a lattice of discrete sites. The 
TASEP is applicable to a broad range of different transport systems, but does not 
consider the fact that in many such systems the availability of resources required 
for the transport is limited. In this paper we extend the TASEP to include 
the effect of a limited number of two different fundamental transport resources: 
the hopping particles, and the "fuel carriers", which provide the energy required 
to drive the system away from equilibrium. As as consequence, the system's 
dynamics are substantially affected: a "limited resources" regime emerges, where 
the current is limited by the rate of refuelling, and the usual coexistence line 
between low and high particle density opens into a broad region on the phase 
plane. Due to the combination of a limited amount of both resources, multiple 
phase transitions arc possible when increasing the exit rate (3 for a fixed entry 
rate a. This is a new feature that can only be obtained by the inclusion of both 
kinds of limited resources. We also show that the fluctuations in particle density 
in the LD and HD phases are unaffected by fluctuations in the number of loaded 
fuel carriers, except by the fact that when these fuel resources become limited, 
the particle hopping rate is severely reduced. 
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1. Introduction 

The totally asymmetric simple exclusion process (TASEP) is one of the fundamental 
models of non-equilibrium statistical mechanics [U El [3] . Essentially a driven diffusion 
model, it has many applications in physics and beyond, including traffic models [3], 
the movement of molecular motors in biological systems [5], and protein synthesis in 
messenger RNA (mRNA) translation [6] [7] . It also belongs to the same universality 
class as some surface growth models [5] . In this paper we study a constrained TASEP 
where finite resources are shared among several lattices. By finite resources we mean 
a constrained number of both particles and "fuel carriers" , whose role is to provide 
the energy needed to the movement of the particles. Molecular motors requiring ATP 
or GTP molecules are an example of such systems occurring in nature. In this paper 
we introduce a new model which includes the finite availability of both resources, in 
contrast to previous works where the effect of having a finite number of a single type 
of resource was studied in isolation [U [TUl HU H21 Q31 H3] ■ As a result, multiple phase 
transitions can occur when varying one of the fundamental parameters of the model 
-the exit rate (3- while keeping the rest of the parameters constant: the system can 
go from a high density regime, to a shock phase, then to a high density phase again, 
visit the shock phase once more, and finally reach a low density phase. This is a novel 
effect that emerges only by combining both limited resources. We use a mean-field 
approach and verify our results by means of Monte Carlo simulations. 

In its most simple form, the TASEP consists of a ID lattice of L sites upon which 
particles can sit, sec figure [TJa). Each site can be occupied only by one particle at 
a time, and particles move from site to site in one direction (say rightward) with 
a hopping rate k. Since particles cannot pass each other, movement requires that 
the downstream site is vacant. A system with open boundaries, as we shall consider 
here, can display rich dynamics with multiple boundary induced phases |15| . Particles 
are allowed to hop onto the lattice with rate a at one end, and off of the lattice 
with rate (3 at the other. For a system with constant internal hopping rate k it is 
possible to solve the steady-state of the system exactly [HI El [TTl [18], whilst the 
full relaxation dynamics have been solved using matrix methods [El [20] . Mean-field 
methods have also been extensively used [TB] , since they are easily tractable and yield 
a good approximation in many cases. There has been much extension of this simple 
model, for example variable hopping rates (site or particle dependent) [21] [22] [23] . 
extended particles which cover more than one site [6j [24] , branching lattices [25] [26] , 
particles which have multiple internal states [27] [28] [29] , as well as a TASEP with a 
constrained reservoir of particles EH [HI [H] ■ 

For the standard TASEP we denote the occupation of the ith lattice site rii = 1 if 
the site contains a particle and ni = otherwise. The system is characterised by the 
steady-state particle current J (the rate at which particles pass any given point on 
the lattice), and the mean site occupancy (often called density) pi = (ni), where (• • •) 
denotes average over realisations of the system (which we assume is ergodic, so this is 
equivalent to a time average). The average density is therefore given by p = L^ 1 £\ pi. 
There arc four possible phases depending on the values of a and /3: the entry limited 
or low density (LD) phase, the exit limited or high density (HD) phase, a maximal 
current (MC) phase where the current depends only on the internal hopping rate, and 
a mixed LD-HD or shock phase (SP). A mean-field approach [16] (which turns out to 
be exact in the L — > oo limit) can be used to calculate J and pi for given a and (3. 
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Figure 1. Schematic diagrams describing the various systems, (a) The TASEP 
with open boundaries in its most simple form. Particles enter at fixed rate a, hop 
at fixed rate k and leave with rate /3. (b) A finite number of fuel carriers can be 
introduced. When a particle hops, fuel from one carrier is used, (c) Complete 
model with finite fuel carrier and finite particles. The entry rate depends on the 
number of particles in the reservoir, and the hopping rate depends on the number 
of loaded fuel carriers. 



The density in the bulk (far from the ends of the system) is given by 
(LD) pld = Qt/fc for a < (3, a < fc/2 , 

(HD) phd = 1 - 0/k for p < a, p < fc/2 , (1) 
(MC) pmc = 1/2 for P, a > fc/2 , 

and the current is always given by J — kp(l — p). The SP, which occurs for 
a = (3 < fc/2, presents an HD region on the right of the lattice and an LD region on the 
left, separated by a boundary which diffuses freely through the lattice. This has often 
been described using a domain wall (DW) theory [3DJ . Due to the free diffusion of the 
DW a time average of the density in SP gives psp = 1/2, but the current depends on 
the density in the LD and HD regions of the lattice, i.e. Jsp = a(l— a/k) = /3(1— 0/k). 

In this paper we consider several TASEPs which share a common finite pool of 
both particles and fuel carriers. The entry rate, which is the same for each TASEP, 
depends on the availability of particles in a common pool (i.e., particles which are not 
involved with any lattice). A model describing several TASEPs sharing a common 
pool of particles has been introduced and thoroughly studied in [9j [TUl [TTj, where 
the authors use the DW theory along with known exact results. In this paper we 
use an alternative recent mean-field (MF) approach that allows us to simplify the 
calculations j!2j . Importantly, we combine this with a model for a finite pool of fuel 
carriers [131 114] which, as noted above, can be viewed as carriers that provide the 
energy which drives the motion, i.e. allowing the particles to hop. Although we 
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consider a fixed number of fuel carriers, we suppose that it takes a finite time to 
"refuel" them with their cargo once it has been used (figure [TJb) shows a schematic 
representation of this model). We show that novel effects arise when both types of 
limited resource are considered, e.g. multiple phase transitions can occur when varying 
the exit rate /3. The outline of the paper is as follows: in section [2] we summarise the 
previous results for the two models separately, before describing in section |3] a mean- 
field model for a system with both a finite pool of particles and fuel carriers which are 
refuelled at a finite rate (figure [ljc)). We then interpret the mean-field model results 
and compare them with results from Monte Carlo simulations. Finally, in section 5] 
we analyse the effect of both limited resources on the fluctuations in the number of 
particles on the lattice. 



2. Finite resources - Review of previous results 

Wc first introduce and describe a system containing multiple TASEPs in which each 
lattice shares the same reservoir of particles; then we present the concept of fuel- 
carriers and the effect of a finite rate of refuelling on the exclusion process dynamics. 



2.1. Finite number of particles 



In this work we analyse a system of M identical lattices of length L. The total number 
of available particles is N , while the number of free particles in the reservoir is N r . 
Since the lattices arc identical and experience the same injection and depletion rates, 
we observe the same phase for each. We can write the total number of particles as 

N = N r + LMp, (2) 

where p is the density on each lattice. The entry rate of the M lattices depends on 
the number of free particles via a saturating function 

~ N r ~ 



a = ao tanh 



ao tanh 



N* 
N 



LMp 



N* 



(3) 



where the constant an gives the entry rate in the limit N r — > oo and is an intrinsic 
property of the lattices yj . Without loss of generality, we fix the normalisation factor 
N* to be LM/2, i.e. the total number of particles used if all the lattices were in the 
MC phase. 

Throughout this paper wc define the different phases according to the values of a 
and j3 and the resulting density p, following |12j . With this choice of nomenclature we 
solve equations for a in terms of an an d N. Since the densities in each phase are the 
same as those in the standard TASEP (equations (JTJ), for a given set of parameters 
(an, ft, N) we find the resulting a which determines the phase; e.g. if a < ft and 
a < k/2 the system will be in the LD phase. By substituting equation © into these 
inequalities, we get a representation of the different phases on the ctQ-ft plane. 

As a consequence of having a finite number of particles, we encounter different 
regimes for small, mid-range and large values of N. We show typical phase diagrams 
for these regimes in figure [2j If N < LM/2, then the HD and MC phases no longer 



| Equation J3]l is consistent with the function used in [9] 1101 1111 . and is relevant, e.g. for the 
application to protein synthesis. 
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Figure 2. Phase diagrams for a TASEP of length L = 500 with a finite number 
of (a) N = 225, (b) N = 475, and (c) N = 600 particles (infinite amount of loaded 
fuel carriers). For small N there are not enough particles to support the HD and 
MC phases, and the SP (coexistence) line opens into a region. For N > LM/2 all 
four phases can be obtained. If N is increased further the HD phase grows at the 
expense of a shrinking SP. 



exist - there are too few particles to support the high density or maximal current 
phases. Instead, there are only two phases: the LD phase and the SP (figure^a)). As 
described in the previous section, the latter occurs when the entry and exit rates are 
equal, i.e., a = j3, and there is coexistence between an LD region and an HD region; 
since a depends on both ao and f3 (through its dependence on the bulk density), the 
line opens into a region on the olq-P phase plane. That is to say, the condition a — j3 
is fulfilled for a certain range of ao [12]- li N = LM/2, then the lattices can support 
an MC phase, and for N > LM/2 there are enough particles for an HD phase to exist 
(ngure[2jb)). As N is increased, the size of the HD phase on the ao-/3 plane increases, 
at the cost of reducing the size of the SP phase (figure EJc)). For N 3> LM the SP 
phase reduces to a line and we recover the original unconstrained TASEP. 

In the unconstrained TASEP within the SP, the LD and HD regions of the system 
are separated by a domain wall (DW) which can diffuse freely across the lattice. 
However, if there is a finite number of particles, in the case of a single lattice the DW 
is pinned to one position [10] (actually the DW fluctuates about its mean position 
like a noisy damped oscillator) . This is because if the DW were to move to the right 
this would increase the number of free particles, increasing the entry rate and therefore 
driving the DW leftwards. Similarly if the DW moves to the left the number of free 
particles decreases, decreasing the entry rate and driving the DW rightwards. The 
opening of the SP line into a region on the cto-(3 plane is possible because a different 
mean position for the DW corresponds to a different value of ao, while keeping a = f3. 
Hence, the system can maintain a = (3 for different values of ao- As detailed in jllj . if 
more than one TASEP is in contact with the same pool of particles, the DW on each 
lattice once again performs a random walk; there is however a pinning of the total 
number of particles on all lattices. 

2.2. Finite reloading time for fuel carriers 

In many of the systems that can be described by a driven lattice gas, the 
energy required for the advancement of the particles is obtained from some kind 
of finite resource. For instance, molecular motors consume ATP molecules, and 
ribosome movement on mRNAs during protein synthesis requires aminoacylated tRNA 
complexes and GTP. The latter case has recently been described in p~3l [14j [31] , and 
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here we briefly review those results, before in the next section combining this with the 
a finite pool of particles model. 

We consider a finite number T of fuel carriers, T of which arc carrying fuel. Every 
time a particle moves, the fuel from one of the loaded carriers is used, and hence, T 
is reduced by one. The hopping rate of particles depends on the availability of loaded 
fuel carriers, and the empty carriers arc refuelled at a rate V. For simplicity, the 
hopping rate is taken to be directly proportional to the loaded fuel carriers, i.e. 

k = aT, (4) 

where a is a constant. Moreover, the rate of refuelling is taken to depend on the 
number T — T of unloaded fuel carriers as 

_ V (f - T) 
b + T-T' 

which has the form of the well known Michaclis-Menten equation in biochemistry. 
The recharging rate is therefore a saturating function of the number of empty carriers 
(T — T) with maximum value Vb and saturation determined by the constant b. Any 
saturating function will give the same qualitative results, but the above formulation 
allows for a straightforward analytical treatment p~3l [14] |j. 

In this section we describe a collection of M identical TASEPs with a finite 
number of fuel carriers (but no constraints on the number of particles); following the 
common mean-field treatment |16j the particle density on the ith site of each lattice 
is given by 

*g = a(l- Pl )-k Pl (l-p 2 ), 

= kpi-i(l -pi) - kpi(l -pi+i), i = 2,...L-l, (5) 
= hp L -x{l-p L )-p PL . 

The inclusion of a finite pool of fuel carriers leads to the additional equation 

dT ( v\ n i m 

~dt = V+T^T ~ 22 (6) 

3=1 

where the sum is over all of the L — l sites which use fuel carriers on each of the M 
lattices. We assume that the particles do not require a fuel carrier to leave the Lth 
site, i.e. the exit rate /3 is constant. In the steady-state we identify the term under 
the sum in ((SJ) as the particle current, and using equation Q we find 

abJ(L-l)M 

k ~ al V -J(L-1)M> [7} 

i.e., the hopping rate is now itself a function of the current. Following |13[ I14j. upon 
solving equations ((5|) in the steady-state we find the four phases as in the original 
TASEP, but now the current and density are given as follows 

(LD) Jld = pld = D(a) for a < a* and a < f3, 

(HD) J HD = J(P), Phd = 1 - V[fi) for p < a* and < a, 

(MC) J MC =a*/2, pmc = 1/2 fora,/3>a*, W 

(SP) J SP = J(a), p SP = l/2 for a = /?<«*, 

§ An alternative model would be to have a hopping rate which is a saturating function of T, and 
then have a constant refuelling rate. This would give qualitatively similar behaviour to the present 
definitions. Our choice is most applicable to protein synthesis, i.e. refuelling due to an enzymatic 
reaction. 
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with L' = L — 1. The behaviour of these functions as a and /? are varied depends on 
the parameters a, b and V /L'M. By considering the steady-state of equation (JB|) and 
noticing that the maximal value that the recharging rate can possibly have is equal to 
Vo, we note that the particle current is limited from above by Vo/L'M. Hence, there 
are substantially two different cases: (i) if Vq/L'M 3> 1, the recharging rate is very fast 
and the particle current is not influenced by it; we recover the results of the original 
TASEP; (ii) if in contrast Vq/L'M <C 1, the recharging of the fuel carriers can limit 
the value of the particle current. Figure [3] shows the current for sets of parameters 
corresponding to each case. In case (ii) (figure EJb)), J(a) shows a sharp change from 
increasing with a, to almost independent of a (though we note that the derivative of 
J (a) remains continuous). The value of J {a) is severely reduced compared to the 
one obtained for case (i) (see figure [3ta)). We refer to the regime where the current 
appears independent of a as a limited resources (LR) regime, since the rate at which 
fuel is used by the particles approaches the rate at which fuel carriers are reloaded. 
Thus the pool of loaded carriers becomes depleted and the hopping rate k reduces. For 
some choices of a and b, the LR regime exists within each of the phases (LD, HD and 
MC). In the LR regime within the LD phase, the sensitivity of the current to changes 
in a or (3 is greatly reduced, whilst the sensitivity of the density is greatly increased. 
In the MC phase, the current is greatly reduced in the LR regime compared to that 
in case (i). For further details see p~3| IT4] . 

The onset of the LR regime depends particularly on the value of the three 
quantities aT, Vq/LM and b/T, the former two controlling at what value of a or /? the 
onset will occur, and the latter controlling the sharpness of the change in behaviour. 
In the rest of this paper we choose b such that there is a sharp onset of LR, and take 
Vo as the control parameter for the fuel carriers, fixing the other parameters. This 
choice not only gives the most interesting dynamics, but it has also been shown to 
be the biological relevant regime in the context of protein synthesis [131 H3] • When 
the onset of LR is sharp we can estimate the value of a or (3 at which this occurs 
by equating the rate of fuel carrier use (approximately aL'M for small a in LD and 
j3L'M for small j3 in HD) and the maximum recharging rate. This gives 



(9) 
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Figure 3. Plots showing the current as a function of a for different parameters. 
In (a) V = 300 s~\ and in (b) V = 20 s -1 . In both cases a = 2 X 10~ 4 , 
T = 5000, and b = 50. Solid lines show the current in the LD phase, Jld = *7(c«), 
and dashed lines the current in the MC phase, Jmc = ct*/2. The dotted line is at 
a*, where there is a transition from LD— >MC. The dot-dashed line in (b) shows 
the value of old as estimated in equation J9]l . 




100 



Figure 4. Steady-state fuelling level T/T as a function of rate Vo for a model 
with no constraints on the number of particles (section l2.2H . Other parameters are 
L = 500, a = 0.5 s" 1 , j3 = 0.1 s -1 , a = 2 X 10~ 4 s" 1 , b = 50, and f = 5000. For 
small Vo the number of loaded fuel carriers is depleted; for large Vo the carriers are 
practically always fully loaded, and we recover the original unconstrained TASEP. 



This value is represented in figure E{b) by a dot-dashed vertical line. As it is shown 
there, the estimation predicts quite accurately the onset of the LR regime. 

3. Constrained reservoir of particles and finite refuelling rate 

A much more realistic model for natural processes such as biological transport has to 
include the finite availability of both particles and fuel carriers. As we show later 
in this section, it is only when combining the two schemes discussed above that 
we can see emerging novel effects, such as multiple phase transitions. Analogous 
to the dependence of the entry rate a on the number of particles N, the steady- 
state proportion of loaded fuel carriers is a saturating function of Vq (see figure 0}. 
Therefore, by regarding N and Vo as control parameters, we can vary the number of 
available particles and loaded carriers respectively; in both cases a saturating function 
of the resource determines the dynamics. 

The quantities of interest are the particle current J, the number of loaded fuel 
carriers T and the particle density p, which is linked to the number of free particles 
by equation ([2]). The expressions for J and p given in equations ([8]) still hold; however 
a is no longer a control parameter, and it can be eliminated using equation (J3]). This 
is the case in all phases except in the SP: since the DW cannot move freely on the 
lattice, the mean density now depends on the size of the LD and HD regions. We 
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calculate the average value pgp using equation ([3]); using the fact that in this phase 
a = P leads to 

PSP^^fl-^tanlW-^V (10) 



LM V N \a 

The hopping rate k (and therefore T) in each phase can be found by substituting the 
appropriate equation for the current J in (|7|). 

We now turn to the problem of finding the boundaries between the different phases 
as functions of ao, 8, N and Vb — i.e., eliminating a. Our aim is to draw the cxq-/3 
phase plane for any given values of Vb and N. As noted in previous sections, due to 
finite particles, the SP line opens into a region, and the HD and MC phases do not 
exist if N < LM/2. Each phase boundary can be written in terms of either ao as a 
function of 8, or vice versa. We now consider the boundaries between each phase in 
turn, consulting figure as an ansatz for the arrangement of the phases. 

(i) MC/LD phase boundary. The MC phase can exist if N > LM/2. If we consider 
starting in the MC phase with large 8 and reducing ao , we will cross into the LD 
phase when a = a*. Using ([3]) this gives an equation for the boundary in the 
ao-/3 plane, where 

a = a* coth [ — ) for 8 > a* and N > LM/2 . 

\N* 2N*J ' 

Note that a* depends on a, b, T, Vb, L and M. Hence the MC/LD boundary is 
a vertical line on the otQ-B phase plane. 

(ii) MC/HD phase boundary. If we now consider starting in the MC phase with large 
ao and reducing B, then we cross into the HD phase when B = a*, i.e., in the 
ao-8 plane this boundary is given by the horizontal line 

8 = a* for a > a* coth ( — | and N > LM/2 . 

(Hi) HD/SP phase boundary. Here we consider moving from the HD to the SP. 
For the system to be in the HD phase we require B < a* , N > LM/2 and 
8 < a. Using equation ^ in the latter inequality gives 8 < ao t&nh(N r /N*); 
using equation ([2]) and the density in HD phase gives the following equation for 
the boundary 

a = /3coth (J^ - ^(1 - V{8))) for 8 < a* and N > LM/2 .(11) 

Hence, the HD/SP phase boundary is a curved line on the ao-8 phase plane. 

(iv) LD/SP phase boundary. Finally we consider starting in the LD phase and moving 
to the SP. In LD we require a < 8 and a < a* . Finding an expression for these 
inequalities poses some difficulty, since by using equation ([3]) and the density in 
the LD phase we obtain 

f N LM 
a = a tanh I — - -j^V{a) 

an equation which cannot be solved analytically to find a as a function of ao- 
Instead we solve this numerically, setting a — 8 (which occurs at the SP) to find 
8 as a function of ao- This gives another curved line on the ao-8 phase plane. 
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Figure 5. Colour on-line. Phase diagrams in the ao-jS plane at different values 
of Vb and N . White lines show the phase boundaries as determined by the mean- 
field model of section [3] Colour maps show the steady-state particle density p 
from Monte Carlo simulations. Parameters used arc V = 20 or 300 s _1 , and 
N = 200, 450, or 600. Other parameters are L = 500, M = 1, a = 2 X 10 -4 s" 1 , 
b = 50, and T = 5000. The dotted line in (f) shows the line of constant ao used 
in figure [8] 



We can then construct the phase plane by plotting the phase boundaries (/3 as a 
function of «o) f° r anv given values of TV and Vo- Unless otherwise stated, throughout 
the rest of this paper we use parameters T = 5000 so that T ;§> LM and hence, it is 
always the refuelling which is the limiting process, and not the total number of fuel 
carriers. This represents realistic scenarios in biological transport processes, such as 
protein synthesis. We set the time scale of the system by choosing a = 2 x 10 -4 s" 1 , 
such that the maximum hopping rate is k = 1 s _1 . A value of b = 50 then gives a 
sharp onset of LR as shown in figure EJb). The phase diagram boundaries calculated 
using the mean- field approach are shown in figure [5] using white lines. 

To test the validity of the mean-field results derived above we perform simulations 
using a continuous time Monte Carlo method [33] . The length of the Monte Carlo time 
step is chosen from an exponential distribution, such that the events occur according 
to a Poisson process, with a single event occurring at each step. Possible events are 
the movement of a particle (cither on to, along, or off of the lattice) or the refuelling of 
a fuel carrier. The event which occurs is chosen stochastically from the set of particles 
which have a vacancy to their right and the set of empty fuel carriers. Particles are 
chosen with a probability such that they move with a rate k, and empty fuel carriers 
arc chosen with a probability such that they are recharged with rate V; after each 
event T is updated accordingly. To remove any transient effects associated with the 
initial condition we disregard the first 5 x 10 6 time steps. Assuming that the system 
is ergodic we average currents and densities over at least a further 4 x 10 7 steps. 
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Figure 6. Colour on-line. Colour maps showing Monte Carlo results for 
the steady-state fuelling level T/T at different values of cto an d P f° r small 
Vfo = 20 s . White lines show phase boundaries as predicted by the mean- 
field model, (a) Small N = 200; (b) mid-range N = 450; (c) large TV = 600. For 
these parameters we have olr = /3lr ~ 0.040 from Other parameters are as 
figure \E\ 



In figure [5] we present a series of phase planes for different values of N and Vq . 
We show the particle density p obtained from Monte Carlo simulations as a colour 
map with the phase boundaries from the mean-field model overlaid (white lines) . Note 
that the mean-field model very closely predicts the boundaries. Here we show data 
for M = 1 lattices, but the plots look the same for M > 1 with appropriately scaled 
parameters. Introducing more lattices does not change the macroscopic behaviour, 
but changes the microscopic behaviour for the SP (see section [3TT]) . 

In figure O we note that for small Vq the phase diagrams look broadly similar 
to the large Vq case, but as we expect from [T3l [14], the phase transitions occur at 
much smaller values of a and /3. We also obtain a limited resources (LR) regime 
within each of the phases. This can seen in figure |H1 which shows colour plots for the 
proportion of charged fuel carriers (T/T) for different values of N, with small V . The 
onset of the LR regime can be clearly seen as T/T decreases dramatically over a small 
range of «o or j3. For large Vq the refuelling is so quick that T/T is constant through 
all phases, i.e. we recover the results for a TASEP with a finite pool of particles, but 
no constraints on the fuel carriers, since these are refuelled almost as soon as they are 
used (data not shown). For clarity, throughout the rest of the paper when we refer 
to LR or limited resources, we specifically mean the regime where the pool of loaded 
fuel carriers has become depleted. 

Crucially, the presence of the LR regime also alters the shape of the phase 
boundaries; a noticeable "kink" can be seen in the LD/SP and HD/SP phase 
boundaries at the point of LR onset. Particularly striking is the shape of the HD/SP 
boundary for the large N small Vq case, and we examine this in detail in section [ 



3.1. The Shock Phase 

As noted above, due to the finite number of particles, the coexistence line present 
in the original TASEP - the shock phase - opens out into a region on the phase 
plane. We examine the behaviour in this phase by locating the domain wall (DW) 
which separates the regions of LD and HD, and examining how the position of this 
is affected by the finite number of fuel carriers and how it changes at different values 
of cto and f3. The introduction of a finite number of particles also gives a change 
from a DW which can wander freely along the lattice, to one where the wandering is 
constrained by the presence of the reservoir. This is most easily explained in the case 
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of M = 1 lattices, where the DW is on average fixed in position. A fluctuation which 
leads to movement of the DW away from its mean position will change the number of 
particles in the reservoir; this in turn changes the entry rate a, acting like a restoring 
force on the DW. In actual fact the DW executes Gaussian fluctuations about its 
mean value, and we discuss fluctuations further in section [4] If multiple lattice are 
introduced (M > 1), then as in [TT] the DW is not pinned, but rather it is the sum of 
the DW position on all lattices which executes Gaussian fluctuations about a mean. 
We focus on the M = 1 case for the rest of this section. 

The relative mean position x E [0, 1] of the pinned DW (where xL gives the site 
at which it is located) can be estimated from the mean-field model by approximating 
the density in the SP as follows 



Since x is always selected such as to maintain the condition a = (3, we can use equation 
(IIOP and the densities from (JSJ to show that 



We also note that the difference between pld and phd decreases as (3 increases, i.e., 
the "height" of the wall decreases. 

As it is the most interesting case, we focus on parameters where the SP has the 
largest area on the phase diagram, namely the mid-range N cases, i.e, figures EJb) 
and (e) . In figure [7J we show plots for the mean position of the DW as a function of 
(3, for large and small values of Vb- Also shown for each case is the density in the HD 
region of the lattice (to the right of the DW) as a function of (3, which will aid in the 
following discussion. 

For large Vo (figures [7] (a) and (c)) there is a monotonic increase in x with j3. A 
larger value of f3 requires that more particles are present in the reservoir in order to 
achieve a = /3. At small j3, x is approximately constant with /3; this is because the 
decrease in density on the HD side of the lattice causes a sufficient release of particles 
to maintain a = j3. Due to the saturating form of the function a(N r ) (equation ([5]) ). 
for larger values of N r a greater increase in N r is required to give the same increase in 
a. Hence for larger values of j3, the change in the HD density as j3 increases no-longer 
releases sufficient particles to keep a = j3; the DW also must move towards the right 
such that there is a steep increase of x with (3. 

In the small Vo case there is an LR regime within the SP which results in an 
interesting dependence of x on j3 (figure. UJTd)); in contrast to the large Vo value case, 
x does not increase monotonically with f3. We can understand this behaviour by again 
considering the density in the HD region of the lattice. We note that for the small Vq 
case, phd changes differently with j3 depending on whether the system is in the LR 
regime or not, and that the maximum in iat/3w 0.04 s _1 corresponds to the onset of 
the LR regime. For (3 < 0.04 s _1 we see from figure [Jjd) that, increasing (3 results in 
a decrease in the density in the HD region - and therefore a release of particles to the 
reservoir and an increase in a. However the decrease of /?hd with (3 is not enough to 
maintain a = (3. The DW must also move rightwards, i.e. there is an initial increase 
of x with j3. After the onset of LR, (3 > 0.04 s _1 , figurc[7Jd) shows that phd decreases 
much more quickly with increasing j3. So now the density on the HD side of the DW 
decreases much more rapidly as (3 increases. The resulting release of particles would 
be too great to maintain a = (3 if the wall did not also move leftwards - x decreases 
again. 



Psp = x Pld + (1 - x)phd- 
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Figure 7. Plots (a) and (b) show the relative mean DW position for systems in 
the SP with Vb = 300 s _1 and Vb = 20 s _1 respectively. In both cases N = 450; 
in (a) a = 0.6 s — 1 , and in (b) a = 0.1125 s _1 . In (b) the onset of the LR 
regime is at /? rj 0.04 s~ x . Plots (c) and (d) show how the density in the HD 
phase phd = 1 — varies with 0, again for Vb = 300 s _1 and Vb = 20 s — 1 

respectively. In the SP this is the density to the right of the DW. In (d), initially 
the density decreases slowly with increasing /?; at the onset of LR the rate of 
density decrease becomes more severe - a change in the behaviour not seen in a 
model with an infinite number of fuel carries. 

Deeper within the LR regime figure [TJd) we have the same situation as before: 
due to the saturating function a(N r ), for large N r we need a greater increase in N r to 
give the same increase in a. The wall has to move rightward as (3 increases in order 
to release enough particles to maintain a = j3. 

3.2. Multiple Phase Transitions 

By combining the effects of both types of limited resources, we obtain a novel phase 
diagram (figure WLX) and figure Etc)). There is an unusual kink shape in the phase 
boundary between the HD and SP regimes, given by equation (fTTj). As can be seen in 
figure Etc) , this is at the point where the system enters the limited resources regime. 
It is possible to draw a vertical line at constant cto through the phase diagram (dashed 
line in figure EJf)), which cuts through the phases HD— >SP— 5-HD— !>SP— >-LD as (3 
increases, i.e., by varying only one parameter we can go from the HD phase through 
a transition to SP, and then a transition back to HD, etc. Figure [5] shows how the 
quantities p, T/T, J, x (where applicable) and a, vary along this line of constant 
«o; we show both Monte Carlo results and the prediction of the mean-field model. 
The mean-field model performs well deep within each phase, but begins to show some 
discrepancy near the phase boundaries; we discuss this further below. 

We label the phases shown in figure [8] with roman numerals I-V and explain each 
in turn. These are that phases which are crossed by the dashed line in figure EJf). 

Phase I At very small values of /3 we have a ^> (3, so the system is in the HD 
phase. Again considering how /?hd varies with /3, from figure EJd) we see that 
for small (3 the slope is small, dpmy/df3 ~ — 1. A decreasing density means an 
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increasing number of free particles, i.e., dN r /d/3 ~ 1; the hyperbolic tangent form 
of equation {3| means that a increases with ft but at a very low rate (da/d/3 <C 1, 
see figure figure [5£d)) [jj 
Phase II We arrive at phase II as follows: in phase I we started with very small 
values of ft such that a>/). By increasing ft, particles are freed and therefore, 
a also increases. However, da/dft <C 1 in phase I, and hence, we eventually reach 
a — ft, and there is a transition to an SP - phase II. Here we have coexistence of 
both LD and HD separated by a DW. 

The current — and therefore the fuel carrier use rate — increases with ft through 
phase I, and initially in phase II (figure Efb)). About halfway through phase II 
(ft « 0.04 s _1 ) the system enters the LR regime (see crosses in figure[5{a)). 
From figure [7£d) we know that in the first half of phase II, the density in the HD 
part of the system decreases slowly with ft — dpnu/dft ~ — 1. The corresponding 
increase in the number of free particles would not be enough to keep a = ft, so 
the DW also moves rightward, i.e. there is an initial increase in x in phase II. 
At the onset of LR the slope of the curve in figure EJd) gets steeper, i.e., 
dp-RT) /dft *C — 1. As ft is increased further (in the second half of phase II) the DW 
must move leftwards again in order to keep a = ft; i.e. after initially increasing 
with ft, x then decreases as LR onsets, as shown in figure [8jc) (see inset). 

Phase III Once the DW reaches the leftmost side of the lattice the system can no 
longer maintain the condition a = ft, so there is a phase transition and we re-enter 
the HD phase. Further increase of ft increases the number of free particles N r , 
however since a is a saturating function of N r , da/dft begins to decrease. That 
is, as ft increases through phase III the slope of a(ft) gets shallower (figure [5Jd)). 

Phase IV If we keep increasing ft, we reach ft = a, and then the system enters the SP 
for a second time. From figure [TJd) we see that deep within the LR regime /?hd 
again varies slowly with ft, and the difference between the LD and HD densities 
is small; therefore in this second SP, changes in the density in the two regions 
of the lattice would not significantly change the number of free particles. Thus 
rapid variation of the DW position is required as ft increases in order to maintain 
a = ft. 

Phase V Once the DW reaches the rightmost edge of the lattice, the system can no 
longer maintain the condition a = ft by moving the wall, and there is a transition 
to the LD phase. 

The above description accounts for the changes in DW position, current and 
density predicted by the mean-field theory, but as we noted previously there is some 
discrepancy with the Monte Carlo results, particularly near the transitions. This is due 
to the fact that our mean-field model assumes that the density is constant throughout 
the lattice, when in fact there is some change near the edges [7]. Also the mean-field 
treatment ignores correlations in the density which occur near the DW. These edge 
effects become less significant as L is increased, and therefore the discrepancy between 
the mean-field and simulation results reduces (data not shown). 

We also note that there is some difficulty in determining the existence and 
position of the domain wall. In figure [5Jc) we define DW position from simulations by 
considering the mean particle density at each lattice site p^, we define the existence 
of a DW if for any pair of adjacent lattice sites i,i + 1 the density cuts through 0.5. 

|| In the model with finite resources a depends on /3, in contrast to the standard TASEP. 
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Figure 8. Colour on-line. Plot showing (a) how the density p (points) and 
fuelling level T/T (crosses), (b) the current J, (c) the mean DW position x, and 
(d) the entry rate a as given by equation (3), vary with /3 for small Vb and large AT 
at a fixed value of ojo = 0-08 s" 1 . Points show Monte Carlo results and lines the 
prediction from the mean-field model. The system passes through several phases 
and we indicate with dotted lines the positions of the boundaries as predicted by 
the mean-field model. We label each phase with roman numerals I-V. The inset 
in plot (c) shows a zoom around the SP phase II. The dashed line in (d) shows 
a = p. 



Then, the position of the DW is given by lattice site i. The difficulty arises in the fact 
that this can also occur near the edges of the system when it is not in the SP. This 
explains why it appears that there are DWs when the system is not in the SP - we 
are actually detecting the decrease in the density at the edge of the system. 
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4. Fluctuations 

In this section we examine fluctuations in the density of particles on the lattice, 
focusing on the case of M = 1 lattice. We obtain power spectra for the fluctuations 
in density by taking the average of the Fourier transform of several different time 
series [32] . These time series are generated by recording the densities at time intervals 
of ~ 125 s. Since our simulation method does not advance time in regular steps, 
the time intervals are not exactly uniform; therefore we use cubic spline fitting to 
obtain a time series at regular intervals. Defining p(t) as the instantaneous density 
at time t, the power spectrum is given by I(oj) = (|FT [p(t) — p] | 2 ), where the angled 
brackets denote average over different simulation runs, and FT[/(t)](w) is the discrete 
Fourier transform of the time series f(t). As before p denotes the time average density. 
From now on, we use the following notation: if we mean the instantaneous value of a 
quantity at time t, then explicit time dependence is indicated; symbols without time 
dependence denote the time average of the quantity. 

4.1. LD and ED Phase 

A good approximation to the power spectrum of the density fluctuations in the LD 
phase for the original TASEP (no particle or fuel carriers constraints) can be found 
via a continuum description, with fluctuations travelling through the lattice with a 
velocity v and an effective diffusion coefficient D [33] At low frequencies (w <C v 2 / D) 
the power spectrum shows oscillations, or dips in power, at unit multiple of 2-kv/L, 
i.e., at frequencies corresponding to the length of time it takes a fluctuation to traverse 
the entire lattice. These oscillations are damped for uj > v 3 /DL, after which I oc cj~ 2 ; 
for large frequencies / oc cj~ 3 / 2 . For the case of finite re-fuelling, but no particle 
constraints (i.e. the N —y oc limit), we obtain a similar power spectrum . As shown in 
figure EJa), we have oscillations with dips at multiples of 2kv/L (with v = k(l — 2p), 
see [34]); at larger values of a, when the system is in the LR regime, the oscillations 
are severely damped, figure [9jb). This is because the hopping rate k is highly reduced 
in the LR regime, meaning that v is reduced, and hence, the damping occurs at 
much lower frequencies. What remains of the dips is still clearly visible in the power 
spectrum at low multiples of 2-kv/L, implying that the density fluctuations are not 
largely affected by the fluctuations in the hopping rate in this regime. The / oc cj~ 2 
for mid-ranged u) and I oc cj~ 2//3 for large ui relationships are maintained even for 
parameters such that the mean number of loaded fuel carriers T is small, conditions 
under which one might have though that fluctuations in T would become important 
(data not shown). We find that actually fluctuations in loaded carrier levels become 
small when the mean value is small 

In summary, just as in the original TASEP, when there are constrained fuel 
carriers the oscillations in the power spectrum of the number of particles on the lattice 
are damped for cj > v 3 /DL. In the original TASEP this is most noticeable at larger 
values of a, which give large values of p, and therefore small v. With constrained fuel 
carriers small v is obtained in the LR regime due to the reduction in the value of k, 
i.e., at much smaller values of a. 

Due to particle- hole symmetry (which is maintained in the present model), our understanding of 
fluctuations in the LD phase can also be applied to the HD phase. 

+ In fact we find that T » ((T(t) - T) 2 > for small T in the LR. The probability of finding T loaded 
fuel carriers appears to be very close to a Poisson distribution in this regime, but further discussion 
of this is beyond the scope of the present work. 
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Figure 9. Plots showing power spectra of the total density from simulations of 
a system with finite fuel carriers with Vo = 20 s — 1 , but no constraints on the 
number of particles (TV — > oo). Plot (a) is for a = 0.01 s — , (b) for a = 0.04 s , 
and (c) a = 0.06 s — 1 ; in the latter two plots the system is in the LR regime, and 
we note that the oscillations have been damped out (see text). The dashed line 
has a slope ui~ 2 , and the dotted lines in each plot are at w = 2nv/L and Attv/L. 
These spectra were obtained by averaging the Fourier transforms of 500 different 
time series. 




Figure 10. Colour on-line. Plots comparing power spectra of density fluctuations 
for systems with a large number of particles (black) and a small number of particles 
(red), (a) Systems not in the LR regime, with Vo = 300 s _1 . Black lines show 
results for a large number of particles: TV = 600 and oq = 0-4 s — 1 . Red lines show 
results for a small number of particles: TV = 220 and ao = 2.83 s — 1 . (b) Systems 
in the LR regime, with Vo = 20 s _1 . Again black lines show results for a large 
number of particles: TV = 600 and ao = 0.07 s _1 . Red lines show results for a 
small number of particles: TV = 220 and ao = 0.65 s _1 . The different values of 
oo used in each case are chosen such that the systems we compare have the same 
mean density and mean hopping rate. 



Density fluctuations have also previously been studied for the TASEP with 
constrained particles [35]. The effect in the LD regime is to suppress the fluctuations. 
This is due to the feedback effect of the particle pool which stabilises the density 
(an increase in density leads to a decrease in entry rate). The effect is reduced for 
higher frequencies, as these correspond to short time scales over which the feedback 
from the particle pool has less influence, i.e. for time scales shorter than the 
time it takes a fluctuation to traverse the entire lattice. Turning to the present 
model with constrained particles and fuel carriers, figure [10] shows power spectra 
for systems with different total numbers of particles, both when there are no limited 
fuel resources ( figure ITDT a)) and when there are (figure fTDT b)). Note that to fairly 
compare fluctuations from two simulations the mean densities and hopping rates (and 
therefore positions of the "dips" in the power spectrum) must be the same; we therefore 
choose different values of ao to give the same value of v in each case. Again we see 
little qualitative difference between the non-LR and LR cases, i.e., fluctuations in the 
number of loaded fuel carriers have little effect on the fluctuations in the density. 
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Jf..2. Shock Phase 

We now consider fluctuations in the density, and in the position of the DW, when the 
system is in the SP. We again consider a single lattice (M = 1), and focus on the case 
of small Vq, and a mid-range value of TV, corresponding to figures [5je) and[6jb), for 
which the SP has the largest area on the ceo-/3 plane. 

To locate the DW we average the occupation of each site over a short time 
window of length t, so as to average over short time microscopic fluctuations in the 
occupation [25] . The value of r is chosen short enough so as to probe the movement 
of the DW on mesoscopic time scales rather than probing only the mean-field density 
profile. Following the approach of [25] we choose for r the smallest value which gives 
exactly one micro-domain wall (/iDW) in the averaged profile. A /iDW is defined as 
any point where the density crosses 0.5 (from above or below) between one lattice site 
and the next. 

We show in figure [TT] the behaviour of the DW at different values of ceo an d 
f3. Figures [TTTa) and (b) show the long time mean density profile for parameters 
for outside and just within the LR regime respectively. Also shown is the mean 
over the mesoscopic time r, i.e. a snapshot of the density profile from which 
the "instantaneous" DW position can be found. We note that in the LR regime 
(figure [TTJb)) there arc larger fluctuations in the r averaged profile than for non-LR, 
and shape of the DW in the long time average profile is wider. The r averaging 
method works well for small values of /3; however for larger /3 when the system is in 
the LR regime, determination of the position of the shock is much more difficult. This 
is because, due to the limited availability of the fuel carriers, the particles move more 
slowly - the density fluctuations in the regions to the left and right of the DW exist 
on time scales similar to that of the movement of the wall. Any value of t which will 
average out the microscopic fluctuations, will also average out the movement of the 
wall. This problem is compounded by the fact that the difference between the mean 
density on either side of the DW decreases with increasing /?. Therefore we can only 
accurately measure the time course of the DW position just inside the LR regime; 
deep within that regime, we can only measure the mean wall position. 

Figures [Tljc)-(d), (e) and (f) show respectively typical time courses, normalised 
histograms of the wall positions, and the correlation function of the time course defined 
as 

_ (Sx(t)6x(t')) t f . 

C[t) - WW) t ' (12) 

where Sx(t) = x{t) — x, and (• ■ ■)( denotes average over time. The correlation 
functions are approximately exponential, and the examples shown have correlation 
time 2.3 x 10 3 s (J3 = 0.007 s" 1 ) and 3.6 x 10 3 s (J3 = 0.049 s" 1 - the LR regime). We 
find that whilst the time scale over which the DW moves is comparable in both cases, 
the width of the distribution is much wider in the LR case. 

For the SP, previous studies [Ml [35] have treated density fluctuations analytically 
by making the approximation that any fluctuations travel quickly along the lattice and 
are absorbed by the DW. That is to say, any fluctuation can be treated as a movement 
of the DW, and so only fluctuations in the rate at which particles move onto the i = 1 
and off of the i = L sites need be considered. The fluctuations can then be described 
using a simple Langevin equation, leading to a power spectra / oc (u> 2 + 7 2 )" 1 , where 
the constant 7 represents the restoring force which localises the DW. 
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Figure 11. Colour on-line. Simulation data showing the behaviour of the domain 
wall in the SP for small Vb and mid-range TV corresponding to figure \5j[e) , for 
ao = 0.1125 s — 1 and various values of 0. Plots (a)-(b) show density profiles, and 
(c)-(d) the DW position time courses for = 0.007 s" 1 (black) and = 0.049 s" 1 
(red). The latter corresponds to a system which is in the LR regime. For the 
density profiles long time averages are shown with heavy lines, and mesoscopic 
time averages (from which we determine the instantaneous wall positions) with 
light lines. The mesoscopic averaging time was r = 500 s. Plot (e) shows DW 
position histograms (normalised) for the same two values of 0, and plot (f) shows 
the time correlation function from 1121 . 




Figure 12. Plot showing power spectra of density fluctuations from simulations of 
systems in the SP with a = 0.1 s" 1 and values of = 0.01, 0.02, 0.03, 0.04, 0.05 
and 0.06 s _1 . The onset of the LR regime is at ~ 0.04 s — , so the top 2 curves 
are for systems with LR. Other parameters are L = 500, Vq = 20 s~~ 1 and TV = 450. 
The dashed line shows the slope u)~ 2 . 
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Figure 13. Plot showing the short time average density profile for a system with 
LR in the SP (heavy line) . The fluctuations arc large and long lived enough that 
they are not averaged out; the DW is almost completely obscured. Also shown is 
the long time average profile (light line) where the DW is visible. Parameters are 
a = 0.1125 s-\ = 0.063 s~\ V = 20 s -1 , and N = 450, and for the short 
time averaged profile r = 500 s. 

Figure [T2l shows that prediction of I oc lj~ 2 for large lj still holds in the system 
with constrained fuel carriers. However, the value of 7 given by the theory does not 
correctly predict the behaviour at small u for the LR regime. As noted above (and 
as can be seen by the increasing values of I{uS) with increasing /3 in figure I12p the 
fluctuations in the density increase significantly in the LR regime. These fluctuations 
are not quickly absorbed by movement of the front, but rather spend considerable 
time in other regions of the lattice; thus a theoretical treatment of fluctuations in a 
system with LR would require consideration of fluctuations in the hopping of particles 
at all sites, and is beyond the scope of the current work. The slow movement of the 
fluctuations along the lattice becomes particularly evident if we look at a snapshot of 
the density profile for a system deep within the LR regime (figure I13|) ; here the front 
seen in the time averaged profile is almost completely obscured by the fluctuations. 

In summary, we find that in the LR regime within the SP density fluctuations 
move more slowly through the lattice. This means that in this regime we cannot use 
the approximation that all fluctuations are absorbed by movements of the DW. 

5. Discussion and conclusion 

In this paper we have introduced and studied a TASEP model which has a constrained 
number of particles, as well as a constrained number of fuel carriers. That is to say, 
there is a finite rate of supply of the energy source which drives the system. 

In a system with only constrained particles O UHl HU [H] , the coexistence or shock 
phase (SP) opens from a line at a = /3 to a region (i.e. a range of ao and j3 values). 
Also, for a very low number of particles (N < LAI/2) the system cannot support an 
HD or MC phase. The introduction of a finite refuelling rate for the carriers leads to 
the existence of a limited resources (LR) regime within the LD, HD, MC and SP. As 
in a system with only constrained fuel carriers [131 114) , the LR regime is reached when 
the rate of fuel carrier use approaches that of the refuelling. 

The main characteristic of the extended model considering both finite particles 
and fuel carriers, is the existence of multiple phase transitions: through increasing 
only the parameter j3, we obtain transitions from an HD phase to a shock phase, then 
back to HD due to the onset of limited resources, then back again to the SP before 
there is a transition to LD. This manifests as a cusp shape on the ao-/3 phase plane. 

Finally we have analysed the fluctuations in the density, and have found them 
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to be broadly in line with those seen in the unconstrained TASEP. For the range of 
parameters studied it appears that fluctuations in the number of loaded fuel carriers 
do not qualitatively change those in the particle density. We do find however that the 
speed at which fluctuations travel decreases in the LR regime, as would be expected 
due to the decrease in the mean hopping rate. This means that the oscillations seen in 
the power spectrum for the LD phase are damped out in the LR regime. The effect of 
a finite number of particles is the same as has been found in previous models [101 135j , 
i.e. the fluctuations are suppressed, particularly for low frequencies. 
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